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Abstract 

Investigation of general-relativistic spherically symmetric steady 
accretion of self-gravitating perfect fluid onto compact objects reveals 
the existence of two weakly accreting regimes. In the first (correspond- 
ing to the test fluid approximation) the mass of the central object is 
much larger than the mass of the accreting fluid; in the second the 
mass of the fluid dominates. The stability of the solutions belonging 
to the first regime has been proved by Moncrief. In this work we report 
the results of a series of numerical studies demonstrating stability of 
massive solutions, i.e., belonging to the second of the aforementioned 
regimes. It is also shown that a formal analogy between "sonic hori- 
zons" in the accretion picture and event horizons in general relativity is 
rather limited. The notion of a "sonic horizon" is only valid in a linear 
regime of small hydrodynamical perturbations. Strong perturbations 
can still escape from beneath the "sonic horizon." 
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1 Introduction 

The model of spherically symmetric self-gravitating accretion in general- 
relativistic hydrodynamics has been studied by Malec [15J and Karkowski et 
al. [9]. It describes a spherically symmetric compact cloud of gas that, due 
to gravity, falls onto a central body (a star or a black hole). 

The underlying mathematical problem consists of solving Einstein equa- 
tions for spherically symmetric space time, where the energy-momentum 
tensor is that of perfect fluid, and where a suitable assumption of quasi- 
stationarity is made — one searches for a solution for which quantities such 
as the density, the pressure, the radial velocity etc., computed at a given 
areal radius, are constant in time. In the accretion case such conditions can 
be satisfied only approximately, as due to the inflow of gas the mass of the 
central object increases and this must be reflected in the change of other 
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quantities (most obviously the space-time metric). Thus, the solutions of 
particular interest are those for which the accretion rate is sufficiently small, 
so that for a long time (as compared to the dynamical scale of the system) 
the central mass remains almost unchanged and the assumption of quasi- 
stationarity can be justified a posteriori. 

It was shown in [9] that for specified boundary conditions such as the 
total mass of the system (the sum of the mass of the central body and 
the mass of the fluid), the areal radius of the accreting cloud and the local 
speed of sound at the outer boundary, there exist two branches of solutions 
that correspond to small accretion rates. Solutions belonging to the first 
branch can be approximated by well known test-fluid solutions (cf. [3] for 
the Newtonian model and [16j for a general-relativistic version) . In this case 
the mass of the central body is much larger than the mass of the accreting 
fluid and the motion occurs practically in the fixed Schwarzschild space-time. 
The other branch consists of solutions for which the mass of the fluid is larger 
than the mass of the central object. In order to obtain these solutions one 
has to solve the complete set of Einstein equations. 

Stability of the solutions in the test-fluid approximation was proved by 
Moncrief [T?J for linearized perturbations. Unfortunately the technique of 
the proof breaks when applied to self-gravitating flows. We have, therefore, 
performed a series of numerical studies of the stability of solutions belonging 
to the massive branch in the full, nonlinear regime. In each examined case 
the simulated flow was stable. 

2 Stationary solutions 

We will consider a spherically symmetric cloud of gas that falls onto a com- 
pact central object. The general, spherically symmetric metric can be de- 
scribed by the line element of the form 

ds 2 = -N 2 dP + Adr 2 + R 2 (d6 2 + sin 2 Odcf 2 ^ , (1) 

where N, A and R are functions of the coordinates t and r. Let us further 
consider a foliation of such a space-time by slices of a constant time t and 
a two dimensional sphere of a constant coordinate radius r embedded in a 
given temporal slice. The trace of the extrinsic curvature of such a sphere 
(twice the mean curvature) reads k = 2d r R/ i^Ry/A^j. 

The accreting gas is described by the energy-momentum tensor of perfect 
fluid 

T" u = {p + p)u»u u +pg» u . (2) 

Here p denotes the pressure, p is the energy density and are the com- 
ponents of the four-velocity of the fluid. The metric tensor is denoted with 
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Let t, r, 9, and appearing in ([T]) be the comoving coordinates. The 
comoving gauge can be obtained by imposing a suitable geometric condi- 
tion expressed in terms of the extrinsic curvature of time slices [15]. Then, 
assuming ([2]), one can show that u r = u e = = 0. 

Let us also introduce a function U = d^R/N. It plays a role of the radial 
velocity in the coordinate system obtained by the transformation (t,r) ^ 
(t' = t,r' = R(t,r)). In these coordinates the components of the four- 
velocity read 

^ = 4 u R = U, u e = u* = 0. 
Another quantity that we introduce is the quasi-local mass 

rRo 



m(R) = m tot ~ 4vr / R pdR' . 
Jr 



Here mtot denotes the total (asymptotic) mass of the considered configura- 
tion, and Roo is the areal radius of the cloud. 
We will also define the local speed of sound a 

2 1 ( , P K 
a = -r [X+ -o 
ft \ n z 



Symbols \ and k denote derivatives 

X : 



dp\ /dp 
dnj e ' \de 

which have to be computed from an equation of state. The quantity n, the 
so-called baryonic density, is defined as a function assuring the conservation 
of the baryonic current V At (nn At ) = 0. Another quantity h = {p + p)/n is the 
specific enthalpy. For barotropic equations of state the above formula can 
be reduced to a 2 = dp /dp. 

By a quasi-stationary solution we will mean a solution for which neither 
the accretion rate defined as m = dfrn nor other hydrodynamic quantities, 
computed at a given areal radius, depend on time. In mathematical terms 
d t 'X = 0, where x = m,U,p, p, . . . For the accretion case such an assump- 
tion can only be satisfied approximately, as we have already mentioned in 
the introductory section. Another assumption on the solution is that of 
asymptotic flatness of the space-time. The accretion cloud of our model is 
compact and has a finite mass. In the simplest version of the model such 
a cloud can be joined to the outside vacuum Schwarzschild space-time by a 
thin transient zone. Another possibility is to consider so-called "quasi-stars" 
[2]. They consist of a spherically symmetric accretion cloud surrounded by 
a huge stellar-type atmosphere providing a reservoir of gas to feed the ac- 
cretion process. 
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Under the aforementioned conditions, the Einstein equations and the 
equations of conservation of the energy-momentum tensor can be written as 



Here B and C are integration constants. These equations are derived for a 
barotropic equation of state, i.e., one for which p = p(p) or, equivalently, 
p = p(n). From now on, we assume the polytropic equation of state p = 



The set of boundary conditions that we choose for above equations con- 
sists of the areal radius of the outer boundary of the cloud Roo, its total 
asymptotic mass m t ot> the value of the local speed of sound at the outer 
boundary a^, and the asymptotic baryonic density n^. These conditions 
do not specify the solution completely. In order to do so one should also 
specify the asymptotic value of the velocity U^. We will, however, adhere 
to a different, more standard approach. Namely, we will search for the tran- 
sonic solutions — those for which in outer regions of the cloud the flow is 
subsonic (the velocity of the fluid is smaller than the local sound velocity) 
and it becomes supersonic in the central parts of the cloud. Such a flow 
passes through the so-called sonic point, where 2U/(kR) = a. Technically, 
this requirement means that one has to find a suitable value of Uco for which 
the flow will be transonic. 

There are two possibilities concerning the treatment of the central ob- 
ject. Suppose we attempt to integrate the equations starting from the outer 
boundary inward, and assume that a transonic solution is obtained. At a 
certain point the optical scalar 9 + = kR/2 + U will vanish and that will 
correspond to the existence of an apparent horizon of a black hole. The 
areal radius of the apparent horizon will by further denoted as Rbh- By the 
mass of the black hole corresponding to a given solution we will understand 
the quasi-local mass enclosed within the apparent horizon tobh = ti(-Rbh)- 
We will also define the mass of the fluid by mfl u id = m to t — "Ibh • 

Alternatively, we can stop the integration earlier and assume that at a 
given point a solid surface of the central object is met. This latter approach 
is used in the model of spherical accretion with radiation. It is assumed that 
the central object does not reflect the inflowing gas, but the corresponding 
energy is emitted in the form of radiation which then propagates through 
the accreting cloud [TP] . 

Fig.Q]shows the dependence of the accretion rate on the ratio m^d/^hot- 
Here all boundary parameters are kept fixed except for the asymptotic bary- 
onic density (the ratio mfl uic j/mtot scales linearly with n^). In most 
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Figure 1: The dependence of the accretion rate rh on the mfl U id- Here 
m tot = 1. 

cases there are two solutions corresponding to a given accretion rate. One, 
for which the mass of the central object is larger than the mass of the ac- 
creting fluid and the other in which the mass of the fluid constitutes most of 
the total mass of the accretion cloud. It can be shown analytically that the 
maximum of the accretion rate is achieved for mfl u id/mtot = 1/3- Solutions 
for which the mass of the accreting fluid is small can be approximated by 
the well known test-fluid solutions. The other branch of massive solutions 
has been recently discovered in [9]. 

3 Stability 

3.1 Numerical codes 

The stability of massive solutions has been investigated numerically using a 
general-relativistic hydro dynamical code. Some details concerning its con- 
struction can be found in [13] . 

The code is based on a modern high resolution shock capturing hydrody- 
namical scheme, where the (possibly discontinuous) solution of the equations 
of hydrodynamics is obtained by considering local Riemann problems that 
appear at interfaces between adjacent grid zones. In some respect, it is sim- 
ilar to a code constructed by Romero et al. [T5]. The implementation of the 
local Riemann solver is that of Donat and Maraquina [7] . The reconstruction 
procedure required to obtain suitable Riemann states uses the "minmod" 
function of Van Leer [12] . 
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The code is formulated in the standard polar gauge with the metric of 
the form 

ds 2 = -a 2 dt 2 + X 2 dR 2 + R 2 (dO 2 + sin 2 6d<j) 2 ) , (4) 



where the lapse a and X are functions of the areal radius R and time t. 

The essential point in the construction of the shock capturing scheme 
is the suitable, conservative formulation of the equations of hydrodynamics. 
Following certain standards established in numerical hydrodynamics com- 
munity we use the Lorentz factor W = au l and the radial velocity v R = 
u R /W (cf. [T]). Then, the equations of hydrodynamics, i.e., (nu^) = 
and VuT^ = can be written as 

d t q_ + -r^d R (aXR 2 F) =aS- {d t InX) q, 



XR 2 

where q denotes a vector of conserved quantities 

q = (D, 5, r) T = (nW, nhW 2 v R , nW(hW - 1) - p 
F is the flux vector 

F = (nWv R , nhW 2 v R v R + p, nW{hW - l)t; fl ) T , 
and I] denotes the source terms 

( ° \ 

s = (nhW 2 v R v R +p)^f~ ( nhw2 -P)^T + ^R 

K -nhW 2 v R ^ - (nhW 2 v R v R +p) y 

Here vr is the covariant radial velocity vr = X 2 v R . 

In the polar gauge, the Einstein equations reduce to just two partial 
differential equations that can be easily integrated by quadratures provided 
that the hydrodynamical source terms are already computed [8j. The first 
of these equations can be easily expressed as 



d R m = 47rir [nhW z - pj = AttR z {D + t), 
where the quasi-local mass m and the metric coefficient X are related by 

The second one provides the radial derivative of the logarithm of the lapse 

d R lna = X 2 + 4ttR (nhW 2 v R v R : + p 
-2 / m , A „T> ( c..R 



' R 2 + 47tR ( Sv +P 



6 



Thus, at each time-step, basing on the known values of hydrodynamic quan- 
tities and metric functions for the previous time, one computes current values 
of hydrodynamic quantities and proceeds to integrate Einstein equations. 

For dynamical calculations we use the perfect gas equation of state 
p = (r — l)ne, where the specific internal density is defined by the rela- 
tion p = n + ne. For smooth and isentropic flows the perfect gas equation of 
state reduces to the polytropic equation of state introduced in the preced- 
ing section. These two equations of state differ in results for discontinuous 
solutions. 

The initial data for the dynamical code consists of five functions of areal 
radius R: v R , m, n and e. These can be obtained by numerical integration 
of equations ((3|) followed by a transition to polar coordinates. The latter is 
achieved using the following relations 



The equations describing quasi-stationary flows were solved using DASPK 
— a solver for a system of differential and algebraic equations developed by 
Petzold, Brown, Hindermarsh, and Li [JJ. In order to find the transonic 
solution we implemented a suitable bisection method taking into account 
that the transonic solution can usually be integrated toward much smaller 
values of R than the other solutions. 

At the inner boundary of the numerical grid the standard outflow con- 
ditions were assumed. This boundary was set outside the apparent horizon 
but always below the sonic point, so that the outflow conditions could be 
easily implemented (it is well known that satisfying such conditions in a sub- 
sonic region causes difficulties and one usually deals with some unwanted 
reflection from the grid boundary). 

The outer boundary was kept fixed using the values obtained from the 
initial solution (inflow boundary). This condition cannot be easily relaxed; 
it is, for instance, also necessary for the stability of the test-fluid accretion. 

The initial solution can be evolved for an arbitrarily long time without 
producing any noticeable changes. Thus, in order to investigate the stability 
of the flow, the initial data was perturbed in the velocity. The perturbation 
was chosen in the form of a bell-shape profile of a sine function restricted to 
one half of its period. 

3.2 Small perturbations 

Fig- presents the evolution of a very small and compact perturbation. Here 
the initial quasi-stationary solution corresponds to the following parameters. 
The exponent in the polytropic equation of state was set to V = 1.4; the 
mass of the whole cloud m to t was normalized to unity; the areal radius of the 
outer boundary of the cloud was equal i?oo = 10 6 ; the asymptotic baryonic 
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Figure 2: Evolution of small perturbations. Here no denotes the density in 
the unperturbed flow. The snapshots show the density contrast (n — no) /no 
in the chronological order. 

density was set to = 1.6 • 10~ 19 , and the asymptotic sound speed was 
such that = 0.1 (here and in all other results we have adopted the 
gravitational units where c = G = 1). For this solution the sonic point and 
the apparent horizon are located at i?* = 0.82 and -Rrh = 0-34 respectively. 
The mass of the fluid constitutes 83% of the total mass. The accretion rate 
is equal fa = 2.6 • 10 -18 . This value is so small that the growth of the central 
mass can be completely neglected during the entire simulation. 
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The subsequent snapshots show the density contrast (n — no) /no, where 
no corresponds to the baryonic density of the unperturbed flow (the unper- 
turbed solution has to be evolved separately in order to subtract numerical 
noise that naturally appears even due to the interpolation of the initial data 
to the numerical grid of the dynamical code). The initial perturbation di- 
vides into two signals: one propagating toward the center and the other that 
propagates outward. The amplitudes of both signals decrease; the first sig- 
nal eventually disappears in the inner boundary, while the second disperses 
until it is no longer distinguishable from the numerical noise. 

3.3 Strong perturbations 

In the case of stronger perturbations the flow is still stable but the overall 
picture can look slightly different. Fig.[3]shows the evolution of the baryonic 
density for such perturbations applied to the same initial solution as before. 
This time the incoming signal forms a shock wave and some part of it is 
reflected from the central regions of the accretion cloud. This is somewhat 
surprising. A small perturbation passing through the sonic point enters a 
zone in which the gas flows onto a center with a velocity larger than the local 
speed of sound, and thus it cannot be reflected outward; in order to make 
this possible such a signal would have to travel through the surrounding 
fluid with a supersonic velocity. Because of that fact a term "sonic horizon" 
has been coined to name the surface enclosing the supersonic region. This 
reasoning breaks in the nonlinear regime of strong perturbations. In this 
case the sonic point can move; new sonic points can be created; others can 
be destroyed. Moreover, the propagation velocity of shock waves that can 
eventually appear is not limited to the local speed of sound. 

The reflection of a strong and wide perturbation is illustrated on Fig. HI 
Here, for clarity reasons, the initial solution was perturbed slightly more 
than in the last case. The solid line shows the graph of Xv R , while the 
dotted line corresponds to —a. The condition for the existence of a sonic 
point expressed in terms used in our numerical code reads jAi^) = a, i.e., 
it corresponds to the intersection of graphs of Xv R and —a. The incoming 
perturbation creates a new sonic point which starts to move toward the 
center of the cloud. At a certain stage the original sonic point is destroyed 
and, after relatively short time, it is replaced with the newly created one. In 
this process, the reflected signal is released and starts to propagate outward. 

Such behavior can be observed for perturbations whose support increases 
to a size larger than the size of the supersonic zone in the steady solution. 
We can, however, show that a strong perturbation which is initially located 
entirely in the supersonic zone can still escape from beneath the "sonic 
horizon." Such case is depicted on Fig.[5j Here, from an initially smooth 
data two shock waves are formed and the one which travels outward can 
eventually pass the original sonic point. 
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Figure 3: Evolution of the perturbed density. The snapshots are placed in 
the chronological order. The profile corresponding to the unperturbed flow 
is depicted with a dotted line. 

These results hint to the limited validity of the formal analogy between 
"sonic horizons" and event horizons known from general relativity (cf . O [6] ) . 
Contrary to the general-relativistic case the notion of the "sonic horizon" is 
restricted to the linear regime of small and compact perturbations. 

Another observation concerning te case illustrated on Fig. [5] is that this 
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Figure 4: Reflection of the incoming signal from the inner parts of the 
accretion cloud. The velocity Xv R is drawn with the solid line, while the 
dotted line depicts values of —a. Intersections of the two graphs correspond 
to the sonic points (see the discussion in text). 

time the perturbation traveling inward is strong but relatively narrow, and 
there is no reflection from the inner boundary. 
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Figure 5: A strong perturbation escapes form beneath the "sonic horizon." 
Here, similarly to Fig.0J the solid line is used for a graph of Xv R and the 
dotted line for —a. 

4 Conclusions 

In this paper we have discussed the stability of the general-relativistic spher- 
ically symmetric accretion onto compact objects. Calculations which take 
into account the self-gravitation of the accreting gas reveal the existence of 
two classes of weakly accreting steady solutions: one which can be well ap- 
proximated by the test-fluid solutions and the other where the mass of the 
fluid is large. Our numerical results suggest that solutions belonging to both 
classes are stable against spherical perturbations (for the test-fluid solutions 
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this can be proved analytically). The stability of solutions in the analo- 
gous Newtonian model has been reported in 1 1 3 (. II 4 j . The Newtonian results 
presented in [T2] show the stability even for non spherical perturbations. 

The advantage of the numerical approach is that we are no longer re- 
stricted to the linear regime of small perturbations. For strong self-gravitating 
perturbations the flow was stable in all examined cases. Investigation of the 
stability against strong perturbations has lead to a remarkable observation 
concerning the sonic points. In this case they are movable and a strong per- 
turbation which propagates against the direction of the flow can easily escape 
from a supersonic region. This shows that the idea of the so-called "sonic 
horizons" makes only sense in the linear regime of small, purely acoustic 
signals. 

Our analysis of the stability is restricted only to the accretion cloud. In 
all numerical simulations we had to ensure that the gas is delivered to the 
system at a small constant rate. To relax this assumption one could explic- 
itly take into account some physical process responsible for such a "feeding" 
of the accretion cloud. One of the astrophysical possibilities is to consider 
the so-called "quasi-stars" described in [2]. Another family of objects to 
which our analysis might be applied consists of Thorne and Zytkow stars 

MM- 

More realistic models of accretion should take into account the radiation 
emitted and transported through the accreting gas. Such models, which take 
into account the self-gravitation of the fluid, were recently analyzed in |10j . 
In essential parts they agree with the results of our, purely hydrodynamical 
approach. 

The model of self-gravitating accretion presented in this paper can be 
also used as a test problem for complex general-relativistic hydrodynamical 
numerical codes. Steady solutions are easy to be obtained; they are stable, 
and the proper treatment of self-gravity is necessary in order to preserve 
them unchanged in the dynamical code. 

Another advantage of this model is that it allows for the investigation 
of the effects caused by self-gravity in a simple, yet not trivial scenario. 
Although the solutions cannot be written in a closed form, many of their 
properties can be obtained analytically (cf. (9J [TT] ) . 
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